clc;clear;
xd = [71.1 72.4 72.4 72.1 71.4 72.0 71.6]';
n = length(xd);
lamda = xd(1:n-1)./xd(2:n);
range = minmax(lamda');
x1 = cumsum(xd);
s = -0.5*(x1(1:n-1)+x1(2:n));
B = [s,ones(n-1,1)];
y = xd(2:n);
u = B\y;
syms x(t);
x = dsolve(diff(x)+u(1)*x == u(2),x(0)==xd(1));
xt = vpa(x,6);
yucel = subs(x,t,[0:n-1]);
yucel = double(yucel);
yuce = [xd(1),diff(yucel)];
epsilon = xd'-yuce;
delta = abs(epsilon./xd');
rho = 1-(1-0.5*u(1))/(1+0.5*u(1))*lamda';